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Abstract 

We  briefly  describe  the  Generalized  Schur  Algorithm  for  the  superfast  solution  of 
positive  definite  Toeplitz  systems  of  equations  and  its  relationship  with  Schur's 
algorithm  and  the  Szego  recursions.  We  then  present  some  experimental  results 
obtained  with  our  FORTRAN  implementation  of  this  superfast  Toeplitz  solver. 
We  will  see  that  the  algorithm  displays  favorable  behavior  in  that  the  growth 
rates  of  the  resulting  residuals  and  errors  are  comparable  with  those  of  the  Szego 
recursions. 


1.    Introduction 


The  problem  of  efficiently  solving  a  linear  system  of  equations  M  x  =  b  where 
M  is  an  Hermitian  positive  definite  Toeplitz  matrix  is  of  considerable  interest  in 
many  areas  of  mathematics  and  engineering.  While  standard  Gauss  or  Choleski 
methods  for  solving  such  a  system  of  order  n  require  0(n  )  arithmetic  operations, 
there  are  several  algorithms  for  obtaining  the  solution  in  0{n~)  operations 
[21,24,5,18,13].  We  refer  to  these  algorithms  as  fast  Toeplitz  solvers.  More 
recently,  algorithms  have  been  presented  for  solving  a  Toeplitz  system  of  equa- 
tions  using    O(?ilog2?))   arithmetic   operations    [11,6,7,15,17,22].     These   superfast 


f  Presented  at  the  International  Conference  on  Linear  Algebra  and  Its  Applications, 
Valencia,  Spain,  September,  1987.  This  research  was  supported  in  part  by  the  National 
Science  Foundation  under  grant  DMS-8704196. 
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Toeplitz  solvers  are  all  based  on  certain  doubling  procedures.  While  they  are  more 
complicated  to  implement  than  the  fast  methods,  they  require  substantially  less 
computational  work  for  sufficiently  large  problems. 

Current  fast  and  superfast  Toeplitz  solvers  are  numerically  unstable,  and 
therefore  unreliable,  when  applied  to  arbitrary  Toeplitz  systems  [8].  However, 
Cybenko  [9,  10]  has  shown  that  the  algorithms  of  Levinson,  Durbin,  and  Trench 
are  numerically  stable  for  the  class  of  positive  definite  Toeplitz  matrices  in  the 
sense  that  accurate  results  can  be  expected  when  the  problem  is  well  conditioned. 
Stability  of  some  superfast  algorithms  can  also  be  expected  for  the  class  of  posi- 
tive definite  matrices  [8].  We  remark  that  some  fast  and  superfast  algorithms  are 
also  potentially  forward  stable  for  other  classes  of  matrices,  for  example,  for 
totally  positive  Toeplitz  matrices. 

The  numerical  reliability  of  fast  and  superfast  Toeplitz  solvers  is  a  topic  of 
ongoing  research.  Since  they  are  relatively  new  on  the  scene  and,  moreover, 
implementations  of  these  algorithms  are  less  tractable  than  those  of  the  fast  algo- 
rithms, numerical  results  for  superfast  Toeplitz  solvers  are  currently  rare. 
Knowledge  of  the  numerical  reliability  is  important  because  superfast  methods 
allow  for  the  solution  of  problems  which  would  require  a  prohibitive  amount  of 
time  for  a  fast  method.  We  are  not  aware  of  any  experimental  results  for  any 
superfast  Toeplitz  solvers  except  for  those  presented  in  [2]. 

Positive  definite  Toeplitz  matrices  play  an  important  role  in  several  aspects 
of  classical  analysis  including  the  classical  theory  of  Szego  polynomials,  the  tri- 
gonometric moment  problem,  and  the  Carathebdory  coefficient  problem  [l,  16]. 
Several  Toeplitz  solvers  are  direct  manifestations  of  algorithms  and  formulas  that 
were  developed  in  the  classical  theory  (see,  e.g.,  [20,  19,  3]).  In  [3],  the  superfast 
Toeplitz  solver  presented  in  [17]  and  [22]  is  described  in  terms  of  a  generalization 
of  Schur's  classical  algorithm. 

In  this  paper  we  give  a  general  description  of  the  Generalized  Schur  Algo- 
rithm and  its  use  in  the  superfast  solution  of  Toeplitz  systems.    We  then  present 
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some  results  of  our  experimental  code  for  the  superfast  solution  of  real  positive 
definite  Toeplitz  matrices.  Refer  to  [3]  for  a  development  of  the  Generalized  Schur 
Algorithm  and  its  relation  with  fast  Toeplitz  solvers,  and  see  [2,4]  for  descriptions 
of  the  implementation  of  the  algorithm  for  complex  and  real  positive  definite  Toe- 
plitz matrices. 


2.    Some  Toeplitz  Solvers 

The  Toeplitz  solvers  we  consider  can  be  divided  into  two  phases:  a  factorization 
phase  and  a  solution  phase.  In  this  section  we  summarize  some  fast  Toeplitz 
solvers  based  on  the  Szego  recursions  and  Schur's  algorithm.  We  then  present  the 
generalization  of  Schur's  algorithm  for  the  superfast  solution  of  a  positive  definite 
Toeplitz  system  of  equations. 

2.1.  The  Szego  Recursions.  The  Hermitian  positive  definite  Toeplitz  matrix 
M=  Mn  +  l  =  [fJ>j-k]?k-=  o  =  M  defines  an  inner  product  on  the  set  of  polynomi- 
als of  degree  at  most  n  by 

<\;,XA->  :  =  (iJ_k 

The  monic  polynomials  \j ,  degXy  =i,  j  =  0,  ■  ■  ■  ,n  that  are  orthogonal  with 
respect  to  this  inner  product  are  the  (monic)  Szego  polynomials  determined  by  M . 
These  polynomials  can  be  constructed  using  roughly  n2  multiplications  and  addi- 
tions by  the  Szego  recurrence  relations. 


-  4  - 
Szego  Recursions: 

for  j  =  0,  1,  .  .  .  ;    n  —  1 

7;+l    :=  -<l,XX;>/<5; 

Xy+iW  :  =  ^  X;  (X)  +  7/+i  X/M 
efy+1:-^(l-bi+il2) 


where  Xy(^)  :  =  ^;  X(l/ ^)  1S  the  polynomial  obtained  by  conjugating  and  revers- 
ing the  coefficients  of  Xj-  The  quantities  7,-  are  called  Schur  parameters;  they 
serve  as  intermediate  quantities  in  the  Szego  recursions  and  are  often  of  physical 
and  mathematical  significance. 

The  Szego  recursions  comprise  the  first  phase  of  several  Toeplitz  solvers.  In 
fact,  they  are  equivalent  with  the  Levinson-Durbin  for  the  solution  of  the  Yule- 
Walker  equations  of  linear  prediction.  Define  R  to  be  the  unit  right  triangular 
matrix  whose  kth.  column  contains  the  coefficients  of  Xk  (P  <  ^  ^n)-  By 
definition  of  the  Szego  polynomials  we  have 

R*  MR  =D  ,     M~l  =  R  D~lR-*, 

where  D  :=  diag[^]0n  and  5k  :=  <XjfcjXjfc>-  Thus,  the  Szego  polynomials  deter- 
mine the  reverse  Choleski  factorization  of  M  .  This  is  the  factorization  that 
underlies  the  solution  phase  of  Levinson's  algorithm  for  the  solution  of  a  Toeplitz 
system  of  equations. 

In  Trench's  algorithm  a  relationship  between  the  entries  of  M_1  is  used  to 
construct  M  from  its  last  column,  which  is  given  by  xn  an^  ^n  •  This  relation- 
ship follows  from  the  Christoffel-Darboux  formula  for  Szego  polynomials  [24,  20]. 
The  matrix  interpretation  of  the  Christoffel-Darboux-Szego  formula  is  the 
Gohberg-Semencul  formula 


where 
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K  AC+1!  =  r,  T/  -  70*  r0 


^o 
^ 


0 

0 

0 

/'(, 

0 

0 

Pn-l 

P<) 

0 

1 

Pn-l 

Pn 

0 

1 

Pn-l 

0 

0 

1 

and  Xn(^)  =  So"/9;'  X-? '  (^n  — 1)-  See,  e-S-»  [12].  Thus,  the  Gohberg-Semencul  for- 
mula provides  a  factorization  of  M  that  is  determined  by  the  last  Szego  polyno- 
mial xn  a]Qd  its  squared  norm  8n.  The  product  M~lb  can  then  be  obtained  in 
O(nlogn)  operations  by  using  Fast  Fourier  Transform  (FFT)  techniques.  This 
provides  an  O(nlogn)  procedure  for  the  second  phase  of  a  Toeplitz  solver. 


2.2.  Schur's  Algorithm.  Another  class  of  Toeplitz  solvere  is  based  on  the  clas- 
sical algorithm  of  Schur  [23].  Let  <f>  be  a  Schur  function,  which  is  a  holomorphic 
mapping  of  the  open  unit  disk  in  the  complex  plane  into  its  closure.  Schur's  algo- 
rithm generates  a  continued  fraction  representation  of  a  given  Schur  function 
<fi  =  (f)0  by  the  successive  application  of  linear  fractional  transformations  (LFTs). 
Given  a  Schur  function  (j>n-\,  define  77i  =  </'n_1(0)  and 

,     /u  1        </>n-l(X)  -In 

<MX)  :=  T" 


X    l-7n*»-iW 

It  is  easy  to  see  that  (f)n  is  also  a  Schur  function.  Furthermore,  if  [7n  |=1  then 
4>n(\)=^n;  in  this  case  <f>  is  a  rational  function  and  the  algorithm  terminates. 
Otherwise,  we  repeat  the  procedure  to  obtain  0n  +  1  from  6n.  In  this  way  we 
obtain  the  Schur  continued  fraction  representation  of  6, 


4>=    Tn(<j)n)   =   tlot2o    ■    ■    ■     otn(<f>n), 


where 


*nW- 


1  +;7nXr 


The  function  Tn  (0)  is  referred  to  as  the  nth  approximant  of  <fi,  and  (pn  is  called  the 
nth  tail  of  </>.  Thus,  Schur's  algorithm  effects  a  sequence  of  elementary  linear  frac- 
tional transformations  on  the  initial  Schur  function  (f),  and  the  Schur  parameters 
7,-  are  intermediate  quantities  that  determine  the  elementary  linear  fractional 
transformations. 

To  implement  Schur's  algorithm  we  write  each  Schur  function  as  the  quotient 

Q?n  2_j  i =0     n   j 

of  formal  power  series,  <pn  =  — —  = : — ,  with  /3n  0  >  0  as  a  partial  nor- 

Pn  2-jj=oPn,j 

malization.  Organizing  Schur's  algorithm  so  that  the  coefficient  pairs  aok^Qk  are 
entered  and  processed  sequentially  leads  to  the  progressive  form  of  Schur  algo- 
rithm, which  is  given  below. 


Schur's  Algorithm: 

for  k  =  1,  2,  3,  •  •  •      while     \"/k  \  <  1 

enter  Qq^-i,  /?o,*-i 

for  j  =1,  2,  3,  •  ■  ■  ,  k-1 


aj,k-j-l 
Pj.k-j 

Ik   =  ak -l,o/ Pk -1,0* 

&.o  =  ft-1.0  (i~  H  I2) 


-%• 


a 


j-l,k-j 
Pj-l.k-3 


If  the  polynomials  a^n),  /?^n^  of  degree  less  than  n  formed  by  the  first  n  terms  of 
Q'0,  /?0  are  input,  then  the  progressive  Schur  algorithm  calculates  a\.n~k\  (3\n~  ' 
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and  7^,  k  =  1,  •  •  ■  ,  n  using  roughly  n    multiplications  and  additions. 

In  contrast  with  the  Szego  recursions,  Schur's  algorithm  provides  a  class  of 
fast  Toeplitz  solvers  based  on  the  Choleski  factorization  of  M,  M  =  L  D  L* ,  with 
L  unit  lower  triangular  (see,  for  example,  [19,  22,  3]).  In  particular,  we  have  the 
following. 

Proposition  2.1.  If  M  =  \^j-k}jlk=o  =  M  is  positive  definite,  then 
Q'q  j  =  —fij+i,  P0  j  =  ft.-,  J  —  0,1,  .  .  .  ,  n—1,  are  the  first  n  coefficients  of 
formal  power  series  aQ,P0  such  that  0O  =  a0/ '/30  is  a  Schur  function. 

Proposition  2.2.    Let  M,  a^n\  /^n)  be  as  in  Proposition  2.1.    Then 

1.  The  Schur  parameters  obtained  by  Schur's  algorithm  applied  to  aft J\ 
PftJ>  are  identical  with  the  Schur  parameters  generated  by  the  Szego  re- 
cursions applied  to  M. 

2.  The  lower  triangular  matrix  L  D  =  [r  ■  k\?  k-o  is  given  by 
rjk  =ft-,y-A-     (*<./)• 

We  remark  that  the  connection  between  Schur  functions  and  Toeplitz  matrices  is 
provided  by  the  solution  of  the  Carat hebdory  coefficient  problem  and  the 
correspondence  between  Carathebdory  functions  and  Schur  functions.  See,  for 
example.  [16,1].  Proposition  2.2  follows  from  by  comparison  of  the  resulting 
Szego  and  Schur  recursions. 

Schur's  algorithm  can  also  be  formulated  in  terms  of  the  LFTs  Tn  .  It  is 
shown  in  [2,  3]  that 


Tn(r)  = 


tn    +  Vn  T 
Vn    +  In  T 


where  f„ ,  X]n  are  polynomials  that  satisfy  the  recurrence  relations 


Vn 


Vn-\     Sn-1 
£n-l     1h-l 


7B 

1 

» 

Co 

= 

0 

1 

(2.1) 


and   |n(X):«X"  fn(l/X).   ?,„  (X)  :  =  XB  5?n(l/X).     We  refer   to  fn    and   t)n.   which 


have  degree  less  than  n  for  n  >  1,  as  the  nth  Schur  polynomials  associated  with 
the  Schur  function  </>0-  The  Schur  polynomials  are  generalizations  of  the  Schur 
parameters  in  the  sense  that  they  determine  compositions  of  the  elementary  linear 
fractional  transformations  that  the  Schur  parameters  determine. 

The  Schur  polynomials  are  also  important  in  that  they  can  be  used  to  obtain 
the  Szego  polynomials. 

Proposition  2.3.  Let  M,  &$  ,  $3  be  as  m  Proposition  2.1,  and  let  f*'7?*  be 
the  kth  Schur  polynomials  determined  by  <y^n\  /3^nK  Then  the  monic  Szego 
polynomials  corresponding  with  M  are  determined  by 


2.3.    The  Generalized  Schur  Algorithm.    Let 

Tn,k   =  *n  +  l°*n+2c    '   '   '    otn  +  k 

be  the  LFT  generated  by  k  steps  of  Schur's  algorithm  applied  to  </>n ,  and  let 
£n  k,r)n  k  be  the  Schur  polynomials  corresponding  with  Tn  k.  Note  that  for  any 
k  >  0  we  have  (f)n  =  Tn  k{(f}n+k)\  ^hat  is,  the  kth.  tail  of  <fin  is  the  (n+£)th  tail 
of  <f>. 

The  Generalized  Schur  Algorithm  is  based  on  a  doubling  procedure  for  gen- 
erating T2n  from  Tn .    It  can  be  described  in  general  terms  as  follows. 


Generalized  Schur  Algorithm. 

To  obtain  ^2n,r]2n  from  a^2n\  (3^2n\ 

0.  Obtain  fn,  ??n  from  a^n\  (3^n\ 

1.  Obtain  Q'i"),  /?(»)  from  ^2n)}  ffin^  ^  ^  using  ^  =  r^), 

2.  Obtain    £n  „,?/„„     from    c4n),  /?(")    as    £n  ,r/n     were    obtained    from 

3.  Obtain  f2n  f^ii    fr°m   £n>   7?n  >  £n,n»V«    US'mg    T2n    =   ^0,n  °  Tn  ,n  • 

The  algorithm  is  started  by  performing  Step  0  directly  by,  for  example,  set- 
ting £:  =  Q'0  o//?o  o  and  7?i  =  1-    More  generally,  we  can  use  Schur's  algorithm  to 

generate  [7;]]  °  for  a  small  value  of  ?i0,  and  obtain  £no,  fjno  from  the  Schur  param- 
eters and  the  recursions  (2.1).  Now  suppose  Step  0  has  been  completed  for  some 
n  >  n0. 

We  have 


and  in  [2,3]  it  is  shown  that 

ft^/n    -0^n    =5„X"    +0(X"+1). 

We  can  therefore  take 

These  equations  allow  us  to  obtain  a[n\  /?<>)  from  a02n),  f3^n\  fn>  *7n  in  Step  1. 

Step  2  is  the  doubling  step.    We  obtain  fnn   and  7?nn  from  ajn)  and  /?•>")  as 
€on  =  Cn  and  7?o,n  =  ?/n  were  obtained  from  q^")  and  0^n\ 


-  10- 

Step  3,  the  composition  of  Tn  and  Tnn,  is  given  by 

S0,2n    ==??0,nSrz,n    ""    £o,n  1h  ,n  > 
7h,2n    =S0,nSn,n    +  r/o,n  ?7n  ,n  • 

Thus,  the  Generalized  Schar  Algorithm  consists  of  various  polynomial  multi- 
plications and  additions  performed  in  a  recursive  manner.  The  multiplication  of 
polynomials  can  be  efficiently  performed  using  standard  FFT  techniques.  This 
results  in  an  0(n\og2n)  algorithm  for  computing  £„  and  r\n ,  where  n  —  2".  The 
Toeplitz  system  of  equations  Mn  +  1  x  =  b  can  then  be  solved  by  forming  xn  from 
£n  >  7?n  by  Proposition  2.3,  and  using  the  Gohberg-Semencul  formula  to  perform 
the  solution  phase  of  the  algorithm  in  O(nlogn)  arithmetic  operations. 

The  recursive  nature  of  the  algorithm  is  easily  handled  in  a  programming 
language  that  allows  self-referencing  subroutines  (e.g.,  PASCAL  or  APL).  In  such 
an  environment  it  is  easy  to  obtain  a  rudimentary  implementation  of  the  algo- 
rithm. However,  in  order  to  obtain  a  FORTRAN  implementation  of  the  algo- 
rithm we  must  consider  the  doubling  procedure  and  resulting  storage  requirements 
more  carefully. 

To  contrast  the  Generalized  Schur  Algorithm  with  Schur's  algorithm,  the 
coefficients  of  ak,  pk  that  are  calculated  by  each  algorithm  are  displayed  in  Figure 
1  (for  n0  =  l).  While  the  classical  Schur  algorithm  generates  aj.n~k\  (5[n~k^,  ^k, 
0<  k  <  n  from  a0  ,  Po  >  ^ie  number  of  coefficients  of  ak,  j3k  calculated  by  the 
Generalized  Schur  Algorithm  is  equal  to  the  largest  power  of  two  that  divides  k, 
unless  this  power  of  two  is  less  than  n0.  Thus,  the  Generalized  Schur  Algorithm 
computes  only  pieces  of  the  Choleski  factorization  of  the  Toeplitz  matrix. 
Nevertheless,  we  can  obtain  every  Schur  parameter;  these  parameters  are  often  of 
physical  and  mathematical  significance. 


FIGURE  1 

The  kth  row  of  squares  represents  the  coefficients  of  ak  and  /3k  that  are  calculated  by  Schur's  algo- 
rithm and  the  Generalized  Schur  Algorithm  (k  =  0,...,n  — 1).  The  filled  squares  represent  the  quanti- 
ties that  are  calculated  during  one  iteration  of  each  algorithm.  In  Schur's  algorithm  the 
coefficients  of  a0,  (30  are  processed  one  at  a  time.  In  the  Generalized  Schur  Algorithm  the 
coefficients  of  a0  ,/?o  are  processed  in  segments,  each  segment  being  twice  the  size  of  the  previous 
segment. 
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3.    Numerical  Experiments 

We  have  a  working  FORTRAN  code  for  the  Generalized  Schur  Algorithm  for 
real  data  based  on  the  procedure  described  in  [4].  Our  code  implements  the  Gen- 
eralized Schur  Algorithm  in  real  arithmetic,  and  uses  dual  split-radix  FFT  algo- 
rithms to  avoid  the  explicit  use  of  the  bit-reversal  permutation.  The  code  uses 
6  n  +  O(l)  real  storage  locations,  and  an  additional  3n/2  real  storage  locations 
for  the  Toeplitz  matrix  and  the  roots  of  unity  (sines  and  cosines)  required  by  the 
FFTs.  The  nth  monic  Szego  polynomial  xn  and  its  squared  norm  8n ,  where 
n  =2U,  are  obtained  using  (8/3)?i  log22n  +  0(n\og2n)  real  multiplications  and 
(16/3)nlog2~n  +  0(?ilog2n)  real  additions.  We  plan  to  make  our  code  available 
to  the  scientific  community  soon. 

In  [4]  it  is  shown  that  xn  j  ^n  can  be  obtained  using  less  than  8  n  log9"n  total 
real  arithmetic  operations.  To  our  knowledge,  this  is  the  smallest  operation  count 
obtained  for  a  superfast  Toeplitz  solver.  Since  the  Szego  recursions  require  more 
than  2?i2  operations,  the  Generalized  Schur  Algorithm  requires  less  computational 
work  for  n  =  2U  >  256.  However,  our  code  does  not  incorporate  the  considera- 
tions discussed  in  [4,  Section  5].  The  operation  count  for  our  code  is  therefore 
somewhat  more  than  the  operation  count  of  [4]  in  the  lower  order  terms. 

When  the  Gaussian  or  Choleski  factorization  is  used  to  solve  the  system  of 
equations  A  x—b,  the  computed  solution  x  satisfies  (A  -\-E)x  =  b  with 
\\E  II  <e<^(n)  ||A  ||.  where  e  is  the  machine  precision  and  <j)(n)  is  a  slowly  growing 
function  of  n.  It  follows  that  the  residual  r  =  b  —  A  x  =E  x  and  the  error 
e  =  A~]E  x  satisfy 

||r||<e^(n)|U||p||,        ||e  ||  <  e#n)  \\A  ||  \\A~l\\  \\x\\- 

Cybenko  [9,  10]  has  shown  that  <f>(n)  is  of  order  n2  for  the  Szego  recursions  (i.e.. 
for  the  Levinson-Durbin  algorithm  for  the  Yule-Walker  equations).  As  the 
Choleski  decomposition  method  also  yields  <j>[n)  of  order  n2.  we  can  conclude  that 
the    accuracy    of   the    Szego    recursions    is    comparable    with    that    of    CholeskVs 
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method.  Accurate  results  can  therefore  be  expected  from  the  Szego  recursions 
when  the  positive  definite  Toeplitz  matrix  is  well  conditioned. 

We  now  present  some  numerical  experiments  to  compare  the  performance  of 
the  Generalized  Schur  Algorithm  with  the  Szego  recursions.  We  are  pleased  to 
report  that  the  results  we  have  obtained  are  so  far  favorable.  In  particular,  the 
Generalized  Schur  Algorithm  displays  roughly  the  same  stability  characteristics 
(e.g.,  the  rates  of  growth  in  computed  residuals  and  errors)  as  the  Szego  recur- 
sions. The  experiments  were  performed  on  the  VAX  11/750  at  Northern  Illinois 
University. 

We  performed  our  experiments  using  the  leading  principal  submatrices  of  two 
classes  of  real  positive  definite  Toeplitz  matrices  Mn+1  =  \v-j-k\jk  =o  =  Mn+i  • 
We  refer  to  these  classes  as  type  A  and  type  B. 

A  (n  ,9)  refers  to  Mn  +  1  where  //,•  =  9J  ,  —  1  <#<1.  These  matrices  have  a  par- 
ticularly special  structure  [14].  The  Schur  parameters  are  given  by  7y  =(— 6)J  ■ 
For  fixed  n  these  matrices  become  ill-conditioned  as  6  tends  to  ±1.  For  fixed  6 
the  condition  numbers  tend  toward  a  limiting  value  as  n  increases.  These  are 
numerically  banded  matrices,  where  the  bandwidth  depends  on  6  and  the  smallest 
positive  machine  number. 

B(n  ,61,62,p)  refers  to  the  autocorrelation  matrix  Mn+1  of  order  n+1  of  a  vec- 
tor z  =  [^-]016n-1  given  by 

n/4 

£j  =r)j  +  E  <x»(iw*)i 

A  =  l 

where  the  0Jk  are  n/4  equispaced  integral  multiples  of  7f/  8n  in  [^1,^2] 
(0  <  9l  <  #2  5s k/  2),  and  t)j  is  randomly  generated  in  [—  p,p].  In  other  words, 
Mn  +  1  =  TT  T,  where  T  is  the  Toeplitz  matrix  formed  by  the  first  n+1  columns 
of  the  circulant  matrix  whose  first  column  is  x.  These  matrices  become  ill- 
conditioned  as  p  tends  to  zero,  and  the  ill-conditioning  accelerates  as  92  —  01 
decreases. 
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We  first  compare  the  accuracy  of  the  Szego  recursions  (Sz)  and  the  General- 
ized Schur  Algorithm  (GS)  for  computing  the  monic  Szego  polynomial  xn  an<^  tne 
Schur  parameters  [Ty]/'.  In  Tables  1  and  2  we  give  the  square  root  of  the  sum  of 
the  squared  errors  (the  2-norm  of  the  error)  as  well  as  the  maximum  absolute 
error  in  the  coefficients  of  xn  an<^  [lj]i  f°r  the  indicated  matrices  of  type  A  and 
type  B.  The  results  of  the  Szego  recursions  in  double-precision  arithmetic  are  used 
as  exact  answers  for  the  error  calculation.  Note  that  the  errors  for  the  two  algo- 
rithms seem  to  grow  at  roughly  the  same  rate  as  n  increases.  In  well-conditioned 
problems,  the  Generalized  Schur  Algorithm  typically  results  in  a  somewhat  larger 
error  than  the  Szego  recursions,  but  usually  this  difference  is  less  than  an  order  of 
magnitude.  On  the  other  hand,  the  Generalized  Schur  Algorithm  occasionally 
yielded  more  accurate  results  in  the  ill-conditioned  problems  of  type  B. 


TABLE 

1 

Ph 

ase  One  E 

,4(4096,0  6 

rrors 

) 

Szego  polynomial  errors 

Schur  parameter  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0  184e-06 

0.219e-06 

0  119e-06 

0.149e-06 

0  802e-07 

0  112e-06 

0  596e-07 

0894e-07 

8 

0315e-06 

0.635e-06 

0.209e-06 

0.328e-06 

0.136e-06 

0  213e-06 

0.782e-07 

0  104e-06 

16 

0353e-06 

0.738e-06 

0.201e-06 

0.328e-06 

0  181e-06 

0.227e-06 

0.782e-07 

0.104e-06 

32 

0.557e-06 

0780e-06 

0.209e-06 

0343e-06 

0.229e-06 

0.227e-06 

0785e-07 

0  104e-06 

64 

0  719e-06 

0.780e-06 

0.179e-06 

0.343e-06 

0297e-06 

0.227e-06 

0.988e-07 

0.104e-06 

128 

0983e-06 

0  780e-06 

0.184e-06 

0.343e-06 

0.371e-06 

0.227e-06 

0  988e-07 

0  104e-06 

256 

0  120e-05 

0  780e-06 

0.238e-06 

0.343e-06 

0436e-06 

0.227e-06 

0988e-07 

0  104  e- 06 

512 

0  149e-05 

0.780e-06 

0.358e-06 

0.343e-06 

0.512e-06 

0.227e-06 

0  988e-07 

0  104e-06 

1024 

0.175e-05 

0.780e-06 

0.477e-06 

0.343e-06 

0579e-06 

0.227e-06 

0988e-07 

0  104e-06 

2048 

0.199e-05 

0780e-06 

0.596e-06 

0.343e-06 

0.639e-06 

0.227e-06 

0.988e-07 

0.104e-06 

4096 

0.225e-05 

0.780e-06 

0.775e-06 

0.343e-06 
A  (4096,0  75 

0.698e-06 
) 

0.227e-06 

0988e-07 

0.104e-06 

Szego  polynomial  errors 

Schur  parameters  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0333e-06 

0.720e-06 

0.238e-06 

0.536e-06 

0.942e-07 

0.321e-06 

0.894e-07 

0298e-06 

8 

0300e-05 

0  112e-04 

0  155e-05 

0.572e-05 

0.780e-06 

0.168e-05 

0.745e-06 

0.991e-06 

16 

0.475e-04 

0.276e-04 

0.187e-04 

0  llle-04 

0473e-05 

0.249e-05 

0.205e-05 

0  991e-06 

32 

0.521e-04 

0.292e-04 

0.197e-04 

0.113e-04 

0.502e-05 

0.252e-05 

0.205e-05 

0991e-06 

64 

0.591e-04 

0.293e-04 

0.197e-04 

0  113e-04 

0614e-05 

0.252e-05 

0205e-05 

0  991e-06 

128 

0  696e-04 

0  293e-04 

0  197e-04 

0  113e-04 

0  736e-05 

0.252e-05 

0205e-05 

0991e-06 

256 

0.835e-04 

0.293e-04 

0  197e-04 

0  113e-04 

0  866e-05 

0.252e-05 

0.205e-05 

0.991e-06 

512 

0986e-04 

0.293e-04 

0.198e-04 

0  113e-04 

0  100e-04 

0.252e-05 

0.205e-05 

0991e-06 

1024 

0  120e-03 

0  293e-04 

0.198e-04 

0.113e-04 

0  118e-04 

0.252e-05 

0  205e-05 

0.991e-06 

2048 

0  140e-03 

0.293e-04 

0  198e-04 

0  113e-04 

0  135e-04 

0.252e-05 

0  205e-05 

0.991e-06 

4096 

0  155e-03 

0.293e-04 

0  199e-04 

0  113e-04 

0.148e-04 

0252e-05 

0  205e-05 

0.991e-06 

A  (4096,0.83) 

Szego  polynomial  errors 

Schur  parameters  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0449e-05 

0  190e-05 

0.286e-05 

0  143e-05 

0  136e-05 

0  786e-06 

0.131e-05 

0  775e-06 

8 

0.369e-04 

0  173e-03 

0.224e-04 

0.918e-04 

0420e-05 

0  147e-04 

0  229e-05 

0  115e-04 

16 

0  180e-02 

0  290e-02 

0.785e-03 

0  114e-02 

0876e-04 

0.775e-04 

0  643e-04 

0  321e-04 

32 

0  231e-01 

0.636e-02 

0598e-02 

0.189e-02 

0474e-03 

0  112e-03 

0  157e-03 

0321e-04 

64 

0309e-01 

0  655e-02 

0  678e-02 

0.190e-02 

0  638e-03 

0  113e-03 

0  I59e-03 

0  321e-04 

128 

0  444e-01 

0.655e-02 

0.678e-02 

0.190e-02 

0894e-03 

0  113e-03 

0  159e-03 

0  321e-04 

256 

0.547e-01 

0.655e-02 

0.679e-02 

0.190e-02 

0  109e-02 

0  113e-03 

0.159e-03 

0  321e-04 

512 

0629e-01 

0.655e-02 

0  680e-02 

0.190e-02 

0  122e-02 

0  113e-03 

0  159e-03 

0  321e-04 

1024 

0  682e-01 

0  655e-02 

0.681e-02 

0  190e-02 

0  131e-02 

0  113e-03 

0  159e-03 

0321e-04 

2048 

0.725e-01 

0.655e-02 

0681e-02 

0.190e-02 

0  140e-02 

0  113e-03 

0.159e-03 

0  321e-04 

4096 

0769e-01 

0  655e-02 

0  682e-02 

0  190e-02 

0  148e-02 

0  113e-03 

0  159e-03 

0  321e-04 

Ph 

TABLE 

ase  One  E 

2 
rrors 

B(4096;20  ',160 

',100) 

Szego  polynomial  errors 

Schur  parameters  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.507e-08 

0314e-07 

0301e-08 

0.310e-07 

0999e-08 

0999e-08 

0.952e-08 

0.952e-08 

8 

0.880e-07 

0  149e-06 

0.659e-07 

0.838e-07 

0357e-07 

0933e-07 

0.242e-07 

0  838e-07 

16 

0.422e-06 

0.850e-06 

0.185e-06 

0.411e-06 

0  134e-06 

0.287e-06 

0888e-07 

0.158e-06 

32 

0.776e-06 

0.878e-06 

0.294e-06 

0.404e-06 

0.208e-06 

0.340e-06 

0888e-07 

0.158e-06 

64 

0.186e-05 

0.918e-06 

0.532e-06 

0387e-06 

0499e-06 

0413e-06 

0.151e-06 

0.158e-06 

128 

0.256e-05 

0  llle-05 

0.841e-06 

0504e-06 

0.587e-06 

0434e-06 

0.151e-06 

0.158e-06 

256 

0.399e-05 

0  138e-05 

0.138e-05 

0.692e-06 

0.683e-06 

0477e-06 

0.151e-06 

0  158e-06 

512 

0  670e-05 

0.181e-05 

0.238e-05 

0.734e-06 

0.781e-06 

0540e-06 

0  151e-06 

0.158e-06 

1024 

0.115e-04 

0  246e-05 

0.441e-05 

0.955e-06 

0.104e-05 

0656e-06 

0.151e-06 

0  158e-06 

2048 

0.205e-04 

0  344e-05 

0.811e-05 

0.966e-06 

0.139e-05 

0919e-06 

0.151e-06 

0  158e-06 

4096 

0.854e-04 

0.270e-04 

0.297e-04 

0.119e-04 

0.547e-04 

0.168e-04 

0.250e-04 

0.616e-05 

5(4096,20  ',160 

*,0.5) 

Szego  polynomial  errors 

Schur  parameters  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.346e-07 

0.587e-07 

0.257e-07 

0.515e-07 

0.328e-07 

0.570e-07 

0.239e-07 

0515e-07 

8 

0.945e-07 

0.474e-06 

0.596e-07 

0.232e-06 

0.625e-07 

0.261e-06 

0.527e-07 

0.232e-06 

16 

0.807e-05 

0.377e-04 

0284e-05 

0  144e-04 

0  116e-05 

0.527e-05 

0.873e-06 

0.430e-05 

32 

0  llle-02 

0.480e-02 

0.406e-03 

0  195e-02 

0.467e-04 

0205e-03 

0.204e-04 

0.937e-04 

64 

0.436e-02 

0.499e-02 

0  172e-02 

0.203e-02 

0.148e-03 

0  312e-03 

0  588e-04 

0.999  e- 04 

128 

0.501e-02 

0524e-02 

0.156e-02 

0.187e-02 

0.185e-03 

0362e-03 

0.588e-04 

0.999e-04 

256 

0.666e-02 

0.848e-02 

0  153e-02 

0.163e-02 

0  238e-03 

0.490e-03 

0.588e-04 

0999e-04 

512 

0.796e-02 

0.182e-01 

0  158e-02 

0.170e-02 

0.280e-03 

0.809e-03 

0588e-04 

0.999e-04 

1024 

0.139e-01 

0316e-01 

0.189e-02 

0.231e-02 

0.428e-03 

0  140e-02 

0588e-04 

0.999e-04 

2048 

0.461e-01 

0.659e-01 

0.303e-02 

0.353e-02 

0.145e-02 

0.311e-02 

0.116e-03 

0  173e-03 

4096 

0.506e-01 

0413e+00 

0  131e-01 
B 

0.104e+00 
(4096,50  ',100 

0.376e-01 

M.o) 

0.221e+00 

0l32e-01 

0  830e-01 

Szego  polynomial  errors 

Schur  parameters  errors 

n 

2-norm 

max 

2-norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.393e-05 

0.409e-05 

0.238e-05 

0.262e-05 

0.169e-05 

0.177e-05 

0.161e-05 

0.173e-05 

8 

0.490e-04 

0.393e-04 

0312e-04 

0277e-04 

0.270e-04 

0.255e-04 

0.206e-04 

0.187e-04 

16 

0.124e-03 

0.307e-03 

0.830e-04 

0.154e-03 

0.585e-04 

0.120e-03 

0.427e-04 

0.631e-04 

32 

0.319e-03 

0.500e-03 

0.140e-03 

0  170e-03 

0.129e-03 

0.224e-03 

0480e-04 

0.757e-04 

64 

0.411e-03 

0.647e-03 

0.120e-03 

0.165e-03 

0  177e-03 

0.333e-03 

0480e-04 

0757e-04 

128 

0.644e-03 

0  890e-03 

0.141e-03 

0.197e-03 

0.245e-03 

0.479e-03 

0.634e-04 

0757e-04 

256 

0.952e-03 

0  129e-02 

0.165e-03 

0.246e-03 

0346e-03 

0.676e-03 

0.634e-04 

0.783e-04 

512 

0.152e-02 

0.231e-02 

0.196e-03 

0.309e-03 

0.530e-03 

0.104e-02 

0.634e-04 

0989e-04 

1024 

0.260e-02 

0  325e-02 

0.231e-03 

0.357e-03 

0904e-03 

0.155e-02 

0.726e-04 

0  1 14e- 03 

2048 

0.303e-02 

0.439e-02 

0.254e-03 

0.396e-03 

0.122e-02 

0.217e-02 

0.726e-04 

0.114e-03 

4096 

0395e-02 

0  641e-02 

0.282e-03 

0.436e-03 

0.173e-02 

0312e-02 

0769e-04 

0.114e-03 

-  14  - 

Our  second  accuracy  experiment  compares  the  Szego  recursions  and  the  < 
eralized  Schur  Algorithm  for  the  solution  of  a  Toeplitz  system  of  equations.    The 
Gohberg-Semencul  formula  is  used  for  the  solution  phase  for  both  algorithms. 

Let   IImUj  and   ||M||2  denote  the  norm  of  the  matrix  M  that  is  induced  by 
the  vector  1-norm  and  vector  2-norm,  respectively.    Then 

max      ||Me;||2<||M||2<V^   \\M\\V  (3.1) 

1  <  ;  <  n 

Let  M  =  Mn ,  and  define  R  :  =  R  D~i/z  so  that  M_1  =  R  R* .    It  is  shown  in  [2] 
that 

\\Re„\\i<\\M-%<n   ||i?e„|||.  (3.2) 


This  gives  bounds  on  the  2-norm  of  M  l  in  terms  of  the  last  column  of  R  (i.e.,  in 
terms  of  the  coefficients  of  the  normalized  Szego  polynomial    Xn-i/vK-i  )• 

After  computing  a  solution  x  of  the  positive  definite  Toeplitz  system,  the 
residual  r  =  b  —  M  x  is  calculated  in  double  precision,  and  the  system  M  e  =  r 
is  solved  (in  0(?i  lg  n )  operations)  to  obtain  an  estimate  e  of  the  error.    Finally, 

the     upper     and     lower     bounds     for     the     two     ratios     ibr  =  — — \-r, — 77-     and 

l|M||||5|| 

•0,  = — are  obtained  using  (3.1)  and  (3.2).    In  order  for  an  algo- 

\\M\\  \\M-l\\\\x\\ 

rithm  to  be  stable,  these  ratios  should  increase  slowly  with  n.  In  Tables  3  and  4 
we  display  the  computed  lower  bounds  on  ||M-1|l2'  the  computed  relative  residu- 
als and  relative  errors,  and  upper  bounds  on  the  ratios  ipr  and  ipe-  These  results 
were  obtained  using  b  =  [y+l].ri=0.  Similar  results  were  obtained  when  the  com- 
ponents of  b  were  set  equal  to  one  and  when  they  were  randomly  generated.  The 
Szego  recursions  were  stopped  at  n  =  2    . 

Note  that  the  sizes  of  the  residuals  and  errors  for  each  algorithm  are  roughly 
the  same,  as  are  the  rates  of  growth  in  these  quantities.  The  upper  bounds  on  ipT 
and  ip    are  slowly  growing  for  both  algorithms.    The  errors  and  residuals  for  the 


TABLE  3 

Toeplitz  Solver  Residuals  and  Errors 


n 

llAC+\  II* 

llrll 

imr 

A(2 
II  £ 

iii 

S,0.6) 

II 

II 

0, 

* 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.449e+01 

0.376e-06 

0.580e-06 

0.353e-06 

0.483e-06 

0.376e-06 

0.580e-06 

0.594e-07 

0.812e-07 

8 

0.493e+01 

0.221e-06 

0.615e-06 

0.253*-06 

0.631e-06 

0.266e-06 

0.740e-06 

0.388e-07 

0.967e-07 

16 

0.495e+01 

0.3S3e-06 

0.444e-00 

0.594e-06 

0.898e-06 

0.545e-06 

0.631e-06 

O.906e-07 

0.137e-OG 

32 

0.495e+01 

0.2S0e-06 

0.457e-06 

0.540e-06 

0.677e-06 

0.449e-06 

0.731e-06 

0.825e-07 

0.103e-06 

64 

0.495e+01 

0.725e-06 

0.117e-06 

0.910e-06 

0.711e-06 

0.125e-05 

0.202e-06 

0.139e-06 

0.109e-0S 

128 

0.495e+01 

0.567e-06 

0.296e-06 

0.120e-05 

0.107e-05 

0.102e-05 

0.529e-06 

0.184e-06 

0.163e-06 

256 

0.495e+01 

0.105e-05 

0.205e-06 

0.171e-05 

0.140e-05 

0.192e-05 

0.37Se-06 

0.261e-06 

0.214e-06 

512 

0.495e+01 

0.150e-05 

0.327e-06 

0.161e-05 

0.697e-06 

0.278e-05 

0.604  e- 06 

0.246e-06 

0.106e-06 

1024 

0.495e+01 

0.2O6e-05 

0.399e-06 

0.222e-05 

0.893e-06 

0.383e-05 

0.743e-06 

0.338e-06 

0.136e-06 

2048 

0.495e+01 

0.288e-05 

0.446e-06 

0.299e-05 

0.934e-06 

0.538e-05 

0.833e-06 

0.457e-06 

0.143e-06 

4096 

0.495e-f01 

0.410e-05 

0.501e-06 

0.418e-05 

0.944e-06 

0.767e-05 

0.937e-06 

0.638e-06 

0.144e-06 

8192 

0.495e+01 

0.595e-05 

0.683e-06 

0.600e-05 

0.109e-05 

0.111e-04 

0.128e-05 

0.916e-06 

0.167e-06 

16384 

0.495e+01 

0.571e-05 

0.590e-05 

0.107  e-04 

0.901e-06 

32768 

0.495e+01 

0.723e-05 

0.738«-05 

0.135e-04 

0.113e-05 

n 

IHC&  \U 

llrll 

JE 

A(215,0.75) 
II  HI 

* 

V- 

Hi 

II 

f 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.279e+02 

0.703e-06 

0.122e-05 

0.434e-06 

0.169e-05 

0.321e-06 

0.557e-06 

0.102e-07 

0.397e-07 

8 

0.647e+02 

0.704e-06 

0.591e-06 

0.278-05 

0.514e-05 

0.316e-06 

0.265e-06 

0.281e-07 

0.520e-07 

16 

0.799e+02 

0.170e-05 

0.163e-05 

0.905e-05 

0.115e-04 

0.967e-06 

0.923e-06 

0.741e-07 

0.942e-07 

32 

0.803e+02 

0.123e-05 

0.425e-05 

0.140e-04 

0.117e-04 

0.959e-06 

0.331e-05 

0.114e-06 

0.954e-07 

64 

0.803e+02 

0.386e-05 

0.285e-05 

0.161e-04 

0.142e-04 

0.402e-05 

0.296e-05 

0.131e-06 

0.116e-06 

128 

0.803e+02 

0.513e-05 

0.293e-05 

0.289e-04 

0.260e-04 

0.681e-05 

0.389e-05 

0.235e-06 

0.212e05 

256 

0.803e+02 

0.922e-05 

0.255e-05 

0.52Oe-O4 

0.504e-04 

0.148e-04 

0.408e-05 

0.424e-06 

0.41  le- 06 

512 

0.803«+02 

0.384e-05 

0.286e-05 

0.267e-04 

0.207e-04 

0.699e-05 

0.521e-05 

0.217e-06 

0.169e-06 

1024 

0.803e+02 

0.610e-05 

0.300e-05 

0.372e-04 

0.315e-04 

0.120e-04 

0.590e-05 

0.303e-06 

0.257e-06 

2048 

0.803e+02 

0.891e-05 

0.288e-05 

0.363e-04 

0.355e-04 

0.183e-04 

0.593e-05 

0.296e-06 

0.289e-06 

4096 

0.803e+02 

0.951e-05 

0.322e-05 

0.317e-04 

0.303e-04 

0.200e-04 

0.679e-05 

0.258e-06 

0.246e-06 

8192 

0.803e+02 

0.692e-05 

0.334e-05 

0.33Se-04 

0.341e-04 

0.148e-04 

0.714e-05 

0.275e-06 

0.278e-06 

16384 

0.803e+02 

0.107e-04 

0.744e-04 

0.231e-04 

0.606e-06 

32768 

0.803e+02 

0.159e-04 

0.738e-04 

0.344e-04 

0.601e-06 

llAC+\  II* 

II- 

II 

A(2 

IU 

1S,0.8) 
II 

V- 

1> 

Tl 

JO 

JL_ 

Hi 

II 

t 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.721e+02 

0.161e-05 

0.265e-05 

0.741e-06 

0.129e-05 

0.448e-06 

0.739e-06 

0.636e-08 

0.110e-07 

8 

0.320e+03 

0.461e-05 

0.546e-05 

0.661e-05 

0.213e-05 

0.932e-06 

O.llOe-05 

0.127e-07 

0.408e-0S 

16 

0.610e+03 

0.462e-05 

0.124e-04 

0.327e-04 

0.826e-05 

0.985e-06 

0.265e-05 

0.329e-07 

0.831e-08 

32 

0.642e+03 

0.764e-05 

0.213e-04 

0.137e-03 

0.190e-04 

0.225e-05 

0.628e-05 

0.131e-06 

0.182e-07 

64 

0.642e+03 

0.112e-04 

0.192e-04 

0.162e-03 

0.400e-04 

0.468e-05 

0.800e-05 

0  155e-06 

0.383e-07 

128 

0.642e+03 

0.427e-05 

0.154e-04 

0.190e-03 

0.123e-03 

0.249e-05 

0.900e-05 

0.181e-06 

0.118e-06 

256 

0.642e+03 

0.126e-04 

0.172e-04 

0.279e-03 

0.291e-03 

0.101e-04 

0.138e-04 

0.267  e- 06 

0.278e-06 

512 

0.642e+03 

0.113e-04 

0.174e-04 

0.264e-03 

0.129e-03 

0.121e-04 

0.186e-04 

0.253e-06 

0.123e-06 

1024 

0.642e+03 

0.876e-05 

0.172e-04 

0.363e-03 

0.264e-03 

0.120e-04 

0.237  e-04 

0.347e-06 

0.253e-06 

2048 

0.642e+03 

0.486e-05 

0.175e-04 

0.392e-03 

0.301e-03 

0.812e-05 

0.292e-04 

0.375e-06 

0.288e-06 

4096 

0.642e+03 

0.817e-05 

0.150e-04 

0.367e-03 

0.295e-03 

0.156e-04 

0.288e-04 

0.351e-06 

0.282e-06 

8192 

0.642e+03 

0.116e-04 

0.178e-04 

0.421e-03 

0.372e-03 

0.241e-04 

0.371e-04 

0.403e-06 

0.356e-06 

16384 

0.642e+03 

0.221e-04 

0.891e-03 

0.482e-04 

0.852e-06 

32768 

0.642e+03 

0.204e-04 

0.906^-03 

0   15S-  n| 

0.867e-06 

TABLE  4 
Toeplitz  Solver  Residuals  and  Errors 


fl(215,20  ',160  \50) 


IK~+\  II: 


imr 


0, 


*t 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.133e+01 

0.824e-07 

0.824e-07 

0.760e-07 

0.760e-07 

0.545e-07 

0.545e-07 

0.538e-07 

0.538e-07 

8 

0.178e+01 

0.211e-06 

0.235^06 

0.136e-06 

0.249e-06 

0.82Se-07 

0.922e-07 

0.701e-07 

0.129e-06 

16 

0.192e+01 

0.372e-06 

0.22  le- 06 

0.112e-06 

0.170e-06 

0.763e-07 

0.453e-07 

0.529e-07 

0.806e-07 

32 

0.200e+01 

0.515e-06 

0.305e-06 

0.468e-06 

0.173e-06 

0.889e-07 

0.526e-07 

0.212e-06 

0.787e-07 

64 

0.203e+01 

0.408e-06 

0.248e-06 

0.198e-06 

0.878e-07 

0.657e-07 

0.399^07 

0.883e-07 

0.391e-07 

128 

0.205e+01 

0.538e-06 

0.389e-06 

0.270e-06 

0.224e-06 

0.840e-07 

0.607e-07 

0.119e-06 

0.988e-07 

256 

0.206e+01 

0.101e-05 

0.736e-06 

0.904^06 

0.650e-06 

0.156e-06 

0.114e-06 

0.396e-06 

0.285e-06 

512 

0.206e+01 

0.395e-05 

0.771e-06 

0.393e-05 

0.677e-06 

0.612e-06 

0.120e-06 

0.172e-05 

0.296e-06 

1024 

0.207e+01 

0.115e-04 

0.117e-05 

0.115e-04 

0.110e-05 

0.175e-05 

0.177e-06 

0.501e-05 

0.482e-06 

2048 

0.208e+01 

0.228e-04 

0.514e-05 

0.229  e-  04 

0.514e-05 

0.339e-05 

0.763e-06 

0.992e-05 

0.223e-0o 

4096 

0.209e+01 

0.501e-04 

0.164e-04 

0.501e-04 

0.164e-04 

0.723  e-05 

0.237e-05 

0.216e-04 

0.708e-0o 

8192 

0.212e+01 

0.102e-03 

0.375e-04 

0.102e-03 

0.375e-04 

0.144e-04 

0.532^05 

0.432e-04 

0.160e-04 

16384 

0.217e+01 

0.196e-03 

0.197e-03 

0.000e+00 

0.265e-04 

0.000e+00 

0.814e-04 

0.000e+00 

32768 

0.614e+01 

0.469e-03 

0.470e-03 

0.000e+00 

0.399e-04 

0.000e+00 

0.520e-04 

0.000e+00 

n  ||A/n-+i  lis 


jn 


B(21S,20  ',160  ',5) 


V\ 


0, 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.158e+01 

0.847e-07 

0.942e-07 

0.964e-07 

0.119e-06 

0.484e-07 

0.539e-07 

0.560e-07 

0.689e-07 

8 

0.382e+01 

0.304e-06 

0.262^06 

0.192e-06 

0.171e-06 

0.732&-07 

0.631e-07 

0.448e-07 

0.398e-07 

16 

0.440e+02 

0.486e-05 

0.314e-05 

0.338e-05 

0.559e-06 

0.137e-06 

0.886e-07 

0.676e-07 

0.112e-07 

32 

0.974e+02 

0.273e-04 

0.272e-04 

0.201e-04 

0.112e-04 

0.696e-07 

0.693e-07 

0.180e-06 

0.100e-06 

64 

0.112e+03 

0.391e-04 

0.285e-04 

0.266e-04 

0.161e-04 

0.810e-07 

0.590e-07 

0.207e-06 

0.125e-06 

128 

0.118e+03 

0.423e-04 

0.363e-O4 

0.299e-04 

0.194e-04 

0.800e-07 

0.687e-07 

0.220e-06 

0.143e-06 

256 

0.121e+03 

0.600e-04 

0.426e-04 

0.471e-04 

0.201e-04 

0.109e-06 

0.777e-07 

0.337e-06 

0.144e-06 

512 

0.123e+03 

0.729e-04 

0.457e-04 

0.662e-04 

0.196e-04 

0.132e-06 

0.827e-07 

0.468e-06 

0.138e-06 

1024 

0.124e+03 

0.914e-04 

0.484e-04 

0.830e-04 

0.189e-04 

0.161e-06 

0.855e-07 

0.582e-06 

0.133e-06 

2048 

0.124e+03 

0.933e-04 

0.456e-04 

0.833e-04 

0.144e-04 

0.160e-06 

0.783e-07 

0.580e-06 

0.101e-06 

4096 

0.126e+03 

0.675e-04 

0.472e-04 

0.500e-04 

0.705e-05 

0.113e-06 

0.787e-07 

0.345e-06 

0.487e-07 

8192 

0.128e+03 

0.607e-04 

0.753e-04 

0.235e-04 

0.530e-04 

0.996e-07 

0.124e-06 

0.160e-06 

0.360e-06 

16384 

0.132e+03 

0.924e-04 

0.703e-04 

0.000e+00 

0.144e-06 

0.000e+00 

0.462e-06 

0.000e+00 

32768 

0.535e+03 

0.922e-03 

0.129e-03 

0.000e+00 

0.850e-06 

0.000e+00 

0.148e-06 

0.000e+00 

IUcA  lls 


Jl*_ 


B(21S,70  ',110  ',1.0) 


*, 


V* 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.276e+03 

0.223e-04 

0.188e-04 

0.513e-05 

0.109e-04 

0.717e-07 

0.606e-07 

0.113e-07 

0.239e-07 

8 

0.771e+04 

0.654e-03 

0.147e-02 

0.273e-03 

0.125e-03 

0.482e-07 

0.109e-06 

0.183e-07 

0.840e-08 

16 

0.819e+04 

0.720e-03 

0.148e-02 

0.332e-03 

0.238e-03 

0.342e-07 

0.704e-07 

0.200e-07 

0.143e-07 

32 

0.872e+04 

0.118e-02 

0.161e-02 

0.757e-03 

0.662e-03 

0.505e-07 

0.687e-07 

0.417e-07 

0.364e-07 

64 

0.926e+04 

0.222e-02 

0.179e-02 

0.141e-02 

0.986e-03 

0.905e-07 

0.728e-07 

0.722e-07 

0.504e-07 

128 

0.939e+04 

0.331e-02 

0.208e-02 

0.270e-02 

0.124e-02 

0.132e-06 

0.831e-07 

0.136e-06 

0.619e-07 

256 

0.953e+04 

0.419e-02 

0.209e-02 

0.394e-02 

0.122e-02 

0.165«-06 

0.821e-07 

0.194e-06 

0.600e-07 

512 

0.959e-t-04 

0.584e-02 

0.196e-02 

0.554e-02 

0.994e-03 

0.226e-06 

0.755e-07 

0.271e-06 

0.483e-07 

1024 

0.963e+04 

0.669e-02 

0.201e-02 

0.645e-02 

0.888e-03 

0.254e-06 

0.759e-07 

0.313e-06 

0.430e-07 

2048 

0.969e+04 

0.749e-02 

0.202e-02 

0.720e-02 

0.826e-03 

0.283e-06 

0.759e-07 

0.348^06 

0.398e-07 

4096 

0.978e+04 

0.825e-02 

0.194e-02 

0.785e-02 

0.739e-03 

0.318e-06 

0.743e-07 

0.375e-06 

0.352e-07 

8192 

0.101e+05 

0.831e-02 

0.178e-02 

0.117e-01 

0.715e-03 

0.323e-06 

0.684c- 07 

0.544^06 

0.336e-07 

16384 

0.103e+05 

0.802e-02 

0.789e-02 

0.333«-06 

0.358e-06 

32768 

0.116e+05 

0.804e-02 

0.277e-01 

0.400e-06 

0.111e-05 

-  15  - 

Generalized  Schur  Algorithm  were  typically  larger  than  those  of  Szego  recursions 
in  the  well-conditioned  problems,  while  the  Generalized  Schur  Algorithm  some- 
times performed  slightly  better  in  the  relatively  ill-conditioned  problems.  These 
differences  were  rarely  more  than  an  order  of  magnitude. 

Thus,  in  the  examples  we  investigated  the  accuracy  of  the  Generalized  Schur 
Algorithm  closely  followed  that  of  Szego  recursions.  Of  course,  while  we  did  not 
encounter  any  problems  where  the  accuracy  of  the  two  algorithms  was  greatly 
different,  more  experimentation  is  needed  to  obtain  a  better  confidence  in  the 
results  of  the  Generalized  Schur  Algorithm. 
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Table  5  shows  average  CPU  times  required  by  the  Generalized  Schur  Algo- 
rithm and  Szego  recursions  for  the  calculation  of  the  nth  Szego  polynomial  and 
the  /i  Schur  parameters.  Note  that  the  Szego  recursions  would  have  used  about  36 
CPU  hours  for  n  =216,  compared  with  about  38  minutes  for  the  Generalized 
Schur  Algorithm.  The  dramatic  relative  efficiency  of  the  Generalized  Schur  Algo- 
rithm for  large  problems  and  the  potential  reliability  of  the  algorithm  as  indicated 
in  the  above  examples  lead  us  to  believe  the  algorithm  will  become  a  popular 
method  for  solving  huge  positive  definite  Toeplitz  systems  of  equations. 


TABLE  5 

Timing 

comparison 

(CPU  seconds) 

n 

GS 

Sz 

256 

3.8 

2.4 

512 

8.1 

8.6 

1024 

17.8 

32.2 

2048 

40.0 

126.0 

4096 

90.7 

505.8 

8192 

200.8 

2030.9 

16384 

441.2 

8068.3 

32768 

1001.6 

32222.3 

65536 

2257.4 

~ 
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